
function[gmm,W,g]=gmm_fun_y_c_avg_eta(candidate_params, candidate_indices, parameters,Yt, vdjQ, t, side, retailer, step)

global  xij_for_investor eta Dmax Imax

%%parmaeters
parameters(candidate_indices)=candidate_params; 

%sort Yt on id and then dealer 
Yt = sortrows(Yt, [1,16,17]); 
    %Note: this is done in the mainfile as well, should not do anything

%number of trades    
Ntrades  = size(Yt,1);

%index for dealers
idxj        = find(Yt(:,17)==1); %index with first entry of each dealer 
dealer_ids  = Yt(idxj,16);
idxd        = find(ismember(Dmax,dealer_ids));

%label parameters
sigma   = Yt(1,15);
xij     = Yt(idxj,14);

if xij_for_investor==1
    xid     = 0;               %set to 0 if psi does not include xid
    xij_shocks = Yt(:,14) ;    %set to 0 if price excludes xij 
    
elseif  xij_for_investor==0
    xid     = xij;             %set to 0 if psi does not include xid
    xij_shocks =  0;           %Yt(:,14);  %set to 0 if price excludes xij 
end 
q       = Yt(idxj,18);

if retailer==0
    munu    = parameters(1);
    sigmanu = parameters(2);
    c       = parameters(3);
elseif retailer==1
    munu    = parameters(1);
    sigmanu = parameters(2);
end 

%distribution functions
f = @(x) normpdf(x,munu,sigmanu);
F = @(x) normcdf(x,munu,sigmanu);

%Back out shocks nu    
theta_long   = Yt(:,5); %midbloomberg price
price        = Yt(:,8); %OTC price 

if  eta==0
    if side==1
        shocks = (price - theta_long) + xij_shocks;
    elseif side==-1
        shocks = (price - theta_long) - xij_shocks;
    end

elseif  eta~=0 && xid==0 && retailer==1 
   
    vdjQ_long = zeros(size(Yt,1),1);
    for t=1:size(Yt,1)  %Note: careful this only works if Yt is sorted!
        idx0            = find(Yt(t,16)==dealer_ids); 
        vdjQ_long(t)    = vdjQ(dealer_ids(idx0));
    end 

    if side==1 %price = eta*value + (1-eta)(theta + nu - xij_shocks) -> nu = (price - eta*value)/(1-eta) - theta + xij_shock
        shocks = (price -eta*vdjQ_long)./(1-eta)- theta_long + xij_shocks;
    elseif side==-1
        shocks = (price -eta*vdjQ_long)./(1-eta) - theta_long - xij_shocks;
    end
end 
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INSTITUIONAL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if retailer==0
    
%use mean theta per dealer (to have a single theta per day)
thetaj = zeros(Imax,1);
for i=1:Imax
    thetaj(i)  = mean(Yt(Yt(:,16)==Dmax(i),5));
end 

%If eta!=0, compute the value of the dealer for the given parameters
if eta~=0 && xid==0 && retailer==0
      estimate=1;
     [~,~,~,~,~,~,~,vdj,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~] ...
         = model_prediction_avg_eta(q, Yt, parameters, side, estimate, 1,[],idxj, 0); 
end 

%initialize so that all outputes have the same number of dealers
psi        =  NaN*ones(Imax,1);

%psi 
if side==1
        psi(idxd,:) =  sigma*log(sum(exp(1/sigma*(xij+q))))-xid;
elseif side==-1
        psi(idxd,:) =  -sigma*log(sum(exp(1/sigma*(xij-q))))+xid;
end 

%into long shape 
psi_long   = zeros(size(Yt,1),1);
theta_long = zeros(size(Yt,1),1);
vdjQ       = zeros(size(Yt,1),1);

for t=1:size(Yt,1)  %Note: careful this only works if Yt is sorted!
    idx0            = find(Yt(t,16)==dealer_ids); 
    psi_long(t)     = psi(dealer_ids(idx0));
    theta_long(t)   = thetaj(dealer_ids(idx0));   
    
    if eta~=0 && xid==0  
    vdjQ(t)         = vdj(dealer_ids(idx0));
    end 
end 

if eta~=0 && xid==0  
    if side==1 %price = eta*value + (1-eta)(theta + nu - xij_shocks) -> nu = (price - eta*value)/(1-eta) - theta + xij_shock
    shocks = (price -eta*vdjQ)/(1-eta)- theta_long + xij_shocks;
    elseif side==-1
    shocks = (price -eta*vdjQ)/(1-eta) - theta_long - xij_shocks;
    end
end 
   
if side==1
       
        b=Inf;
        a      = (psi_long + c - eta *(vdjQ + xij_shocks))./(1-eta) -theta_long;
           
        Mills   = f(a)./(1 - F(a));  
        ExpfunT = munu + sigmanu.*Mills ;  
        VarfunT = (sigmanu.^2 .*(1 + (a .*f(a))./(1 - F(a)) - (Mills).^2));

        massOTC  = F(b) - F(a);
        massE    = F(a);
        
elseif side==-1 
        
        b = (psi_long - c + eta *(vdjQ + xij_shocks))./(1-eta) -theta_long ; %NOT UPDATED!!
        a= -Inf;
        
        Mills   = - f(b)./F(b);  
        ExpfunT = munu + sigmanu.*Mills ; 
        VarfunT = (sigmanu.^2 .*(1 + (-b .*f(b))./(F(b)) - (Mills).^2));

        massOTC  = F(b) - F(a);
        massE    = 1-F(b);
end   
    
        probOTC  = massOTC./(massOTC+massE); %probability to choose OTC conditional con trading      
        shareOTC = size(Yt(Yt(:,9)==1,:),1)./Ntrades;
       
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% RETAILER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

elseif retailer==1
    
    if side==1
        
        ExpfunT = munu;
        VarfunT = sigmanu.^2;
        
    elseif side==-1 
        
        ExpfunT = munu;
        VarfunT = sigmanu.^2;
    end  
    
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% MOMENTS & GMM function 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

moment1 = (ExpfunT - shocks);
moment2 = (VarfunT -(shocks -ExpfunT).^2);

if retailer==0
    
    idxE = find(Yt(:,9)==2);
    moment1(idxE)=0;
    moment2(idxE)=0;
    
    moment3 = probOTC - shareOTC;
    moment = [moment1, moment2, moment3];
    
elseif retailer==1
    moment = [moment1, moment2];
end 

g =  1/Ntrades *sum(moment);

if step==1 
    W=eye(size(g,2));
elseif step==2 
    W=inv(1/(Ntrades) * (moment') * moment);
end 

gmm0 = g* W* g';

if isnan(gmm0)==0 
        gmm = gmm0;
   else 
        gmm = 10^100;
end 
    
  






